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Polyethylene under tensile load: strain energy storage and breaking of linear and 
knotted alkanes probed by first-principles molecular dynamics calculations 

A. Marco SaittaQ and Michael L. Klein 
Center for Molecular Modeling, Dept. of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6202 

(February 1, 2008) 

The mechanical resistance of a polyethylene strand subject to tension and the way its properties are 
affected by the presence of a knot is studied using first-principles molecular dynamics calculations. 
The distribution of strain energy for the knotted chains has a well-defined shape that is very different 
from the one found in the linear case. The presence of a knot significantly weakens the chain in 
which it is tied. Chain rupture invariably occurs just outside the entrance to the knot, as is the case 
for a macroscopic rope. 
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I. INTRODUCTION 

Structural isomerism is a well-studied phenomenon in 
chemistry, but it is also possible to construct isomers of 
molecules that differ not in their connectivity, but in 
their topology ||-||]. These topological isomers exhibit 
different optical rotary power. Examples of topologically 
complicated structures are knots, for which the simplest 
is the "trefoil". A trefoil is constructed from one chain 
that forms a single, self-threaded loop Technically, 
the ends of the chain must be connected for this struc- 
ture to be a knot DNA fragments, for example, have 
been observed to form knots similar to the trefoil as well 
as more complicated ones [Q . Interpenetrating entangled 
polymers also form knots. Statistical arguments based on 
self-avoiding random walks, Hamilton cycles, and so on, 
have demonstrated conclusively that the chain length is 
the crucial factor governing the probability of knot forma- 
tion 1^,^. As the chain length increases, this probability 
goes to one exponentially 

The self- entanglement of polymer chains has profound 
effects on the physical properties of the bulk material, 
such as its crystallization behavior, elasticity, and rheol- 
ogy. For example, the transfer of a knot from a polymer 
in the melt state to that in the crystalline/glassy state 
has been studied . Crystallization of polymers is in- 
hibited by the presence of knots and the freezing tends to 
remove them. On the other hand, quenching the melt to 
a glass tightens and stabilizes existing knots. The struc- 
ture and stability of knots in polymers has been a subject 
of extensive study ]l^-|l5t. We recently reported a first- 
principles molecular dynamics study of a trefoil knot in 
a polyethylene strand, which revealed a surprising simi- 
larity in behavior between macroscopic and microscopic 
knotted strands under tensile stress |]l6|. In the present 
work, we focus on the storage of strain energy in a poly- 
mer. In particular, we will describe how strain distribu- 
tions differ between linear and knotted polymer chains. 
Further, this work demonstrates that an efficient combi- 
nation of classical and density-functional based first- 
principles 111] molecular dynamics (MD) can be a very 



useful tool in the description of physical as well as chem- 
ical properties of polymers with and without topological 
defects. 

The essential details of the calculations are provided 
in the following section. The results for the linear 
polyethylene-like (alkane) molecules are reported in Sec- 
tion III, while classical and first-principles simulations 
on the knotted strands are presented in Section |l^. We 
discuss our findings in Section 

II. COMPUTATIONAL METHOD 

The first principles calculations were performed using 
the Car-ParrincUo (CP) method which combines 

MD and density- functional theory (DFT). DFT is known 
to provide an accurate and reliable description of carbon- 
and hydrocarbon-based systems ||l9| . We adopted Becke 
and Lee- Yang-Parr (BLYP) gradient corrections to 
describe the exchange and correlation functionals, re- 
spectively. The electron-ion interaction was described 
by pseudopotentials in the Martins- TrouUier form. 
The plane-wave basis set had a cutoff of 60 Ry, with a 
F-point integration in the Brillouin zone. Electron spin 
was explicitly taken into account [p2| , to obtain more re- 
liable values of the dissociation energies. We described 
the isolated molecules with periodic boundary conditions 
and supercells that were large enough to avoid any signif- 
icant interaction with periodic images. The same, large, 
tetragonal supercell (12 Ax 12 Ax 20 A) was used in both 
the linear and knotted cases. While these dimensions 
may be oversized for the linear molecules, this approach 
afforded us consistent and comparable results. 

Classical MD simulations were carried out according to 
the united atom (UA) scheme , in which every methy- 
lene {—CH2—) or methyl {—CH^) group is considered a 
"pseudoatom" . The potential of interaction included har- 
monic bond-stretching, harmonic bond-bending, a tor- 
sional field, and a Lennard- Jones interaction between 
non-bonded atoms. We fitted the stretching and bending 
interaction constants of the model to ad hoc CP calcula- 



tions on small alkanes. 

In both classical and first-principles MD simulations, 
external tensile stress was mimicked by constraining the 
chain end-to-end distance L to a fixed value. The system 
was evolved dynamically at room temperature before be- 
ing annealed in order to find the (constrained) minimum. 

III. LINEAR ALKANES 

All the results reported in this section refer to CPMD 
simulations. As a first step, we studied a decane molecule 
(Cioi?22), which can be taken as representative of longer 
linear alkanes and, by extrapolation, of polyethylene. 
The equilibrium structural parameters, reported in Ta- 
ble 1, are in excellent agreement (within 1%) with other 
DFT results and with experimental data |2J,|2^ for hy- 
drocarbons and/or crystalline polyethylene (PE). 

The molecule was stretched in a series of steps. Ad- 
ditional external load was applied to the system by suc- 
cessively increasing L. To do so, the nuclear coordinates 
were all scaled uniformly along the molecular z axis. In 
the first stages of tensile loading, the geometry of the 
system was optimized by minimizing the forces acting 
on atoms, but no dynamics was used. Quite unexpect- 
edly, the strain did not distribute uniformly along the 
chain, but tended to concentrate equally on the two ex- 
tremal bonds, where the force originating from the con- 
straint was actually applied. This behavior can be clearly 
observed by comparing the terminal CCC bond angles 
and C-C bond lengths with the average values calcu- 
lated for atoms in the interior of the chain. These re- 
sults are reported in panels a and b of Fig. |l| as a func- 
tion of the distance L. For deformations AL larger than 
10% of the equilibrium length Lq, we followed the pro- 
cedure, described in the previous section, of letting the 
molecule evolve dynamically at room temperature, and 
then cooling it down with a simulated annealing tech- 
nique. The global constrained minimum of the system, 
after the CPMD evolution, was found by geometry opti- 
mization of the structure. The linear n-decane does not 
break for distortions up to AL ~ 18.0%. 
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TABLE I. Equilibrium structural parameters of n-decane 
as determined in the present work (DFT-BLYP) or using a 
Becke-Perdew (DFT-BP) functional for the exchange and cor- 
relation (23), compared to theoretical (DFT-BP) and experi- 
mental results for crystalline polyethylene (PE) (23,24). 



The distribution of strain energy of the system was 
estimated by adopting the UA model to analyze 
the configurations provided by the first-principles cal- 
culations. In the present case, contributions to the to- 
tal strain energy originating from the torsional potential 
and the long-range non-bonded interactions are negligi- 
ble. Bond stretching and bending largely dominate in 
the storage of tensile load in linear molecules. The sys- 
tem can sustain about 16.2 kcal/mol per C-C bond before 
rupture (Fig. ||, panel c). 



140 




AL(%) 

FIG. 1. Analysis of the strain distribution along the 
n-decane molecule as function of the linear deformation. In 
panel a shown is the bond angle (degrees) of the terminal 
carbons (diamonds, solid line), and the average other C-C-C 
bond angles (circles, dotted line). In panel b we report the 
C-C bond-length (A) distribution accordingly. In panel c the 
average strain energy (kcal/mol) per C-C bond is displayed. 

The analysis of the strain energy distribution along the 
chain for L = 1.18 Lq, shown in Fig. ^, confirmed the re- 
sults previously discussed, i.e., the stress was not uniform 
along the molecule. Instead, it was mostly concentrated 
on the two extremal bonds, which store 20-25 kcal/mol of 
the energy associated with the distortions from the ideal 
geometry. All the other "central" bonds sustain a lower 
stress, whose average energy is about 12 kcal/mol. 
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FIG. 2. Strain energy distribution along the C-C bonds 
of the n-decane molecule for a deformation AL = 18%, esti- 
mated by making use of the UA model (see text) . The stress is 
seen to concentrate most on the extremal bonds, while being 
almost constant along the others. 




FIG. 3. A strained (L = 1.18 Lo) n-decane molecule before 
(left) and after (right) bond breakage. Carbons are displayed 
in dark gray and hydrogens in white. The sp^ — > sp^ change 
of the hybridization status of the atoms involved is detectable 
from the coplanarity of the two carbons with their bonds. 

The qualitative features of the strain energy distribu- 
tion were essentially the same for any value of length 
deformation, AL. The bonds closest to the two extremal 
ones typically were about 6% energetically less strained 
than the average of the central bonds. Thus, from the 
third atom on, the system can be considered to be "bulk" , 
and the strain practically constant. Additional calcu- 
lations performed on a linear C2i-ff44 molecule indicate 
that the shape of the strain energy distribution reported 
in Fig. II for C10H22 is unaffected by the size of the sys- 
tem. The amount of strain energy on the extremal and 



the "bulk" bonds are within the errors of the values ob- 
tained for the shorter alkane. 

As foreshadowed by the classical strain energy dis- 
tribution, the CPMD simulation results demonstrated 
that, for larger elongations, the alkane molecule actu- 
ally breaks at one of the extremal bonds. As shown in 
Fig. H, the CgHig- radical contracts very quickly after 
the break. While large deviations from equilibrium bond 
lengths were found throughout the simulations, the true 
marker for bond dissociation was the change in hybridiza- 
tion (sp^ —^ sp^) for the previously bonded carbon atoms. 
Indeed, as shown in the right side of Fig. ||, the three 
bonds about those carbon atoms become coplanar after 
the break. The length of the bond formed between the 
last carbon atom of the CgHig- radical and its closest 
neighbor decreases as well, oscillating around the typical 
C(sp3)-like value, ca. 1.49 A. The experimental enthalpy 
of dissociation of a methyl group from smaller n-alkanes 
is ca. 86 kcal/mol. Our theoretical dissociation energy of 
83 kcal/mol, calculated as the energy difference between 
the optimized C10H22 molecule and the optimized rad- 
icals, gives a very good quantitative description of the 
dissociation phenomenon. 
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FIG. 4. The n-undecane molecule after bond breakage. As 
displayed, the constraint on the distance L is applied between 
Ci and Cio. The resulting sp^ hybridization of C9 and Cio is 
clearly observable. 

The finding that a linear alkane subject to tensile 
stress breaks where the external force is applied was 



checked by performing analogous simulations with a n- 
undecane molecule (CiiiJ24). In this case, however, the 
constrained distance L was not the end-to-end length of 
the chain, but still the Ci-Cio distance. Indeed, this sys- 
tem behaved very similarly. As was observed for the n- 
decane molecule, the chain once again ruptured where 
the constraint was applied. In contrast to the previ- 
ous case, however, the simulation consistently formed an 
ethyl radical (Ciii?24 ^ CgiJig- + C2Hq-) (see Fig. |). 
The present DFT value for the dissociation energy is 
81 kcal/mol, which is again in excellent agreement with 
the experimental dissociation enthalpy of an ethyl group 
which is about 82 kcal/mol. We note here that similar 
calculations that do not take into account spin variables 
give the same qualitative results, but the corresponding 
dissociation energies (102 and 94 kcal/mol, respectively) 
are overestimated by 10-15 %. 



rectly involved in the formation of the knot and hence 
they could be removed (Fig. ||, panels b and c). 

Longer simulations were performed when the size of the 
molecule was reduced to iV = 50 pseudoatoms. First, the 
system was kept at room temperature for at least 200 ps 
before cooling it down, with a slow annealing (about 50 
ps), to r = -ftT. The strain-energy distribution of each 
chain at the constrained (knotted) minimum energy con- 
figurations was systematically determined for each value 
of the external load. The results obtained from the CP 
calculations on the linear hydrocarbons suggest that such 
systems can sustain a strain-energy per bond in the range 
of 15-30 kcal/mol before rupture. It is thus noteworthy 
that the per bond strain energy for iV = 50 chain is negli- 
gible on this "relevant" energy scale (Fig. even though 
the shape of the trefoil is very well defined (Fig. ^, panel 
c). 



IV. KNOTTED MOLECULES 
A. Classical MD calculations 

The study of a knotted polymer requires a non-trivial 
procedure to set up the starting nuclear coordinates for a 
trefoil. Ideally, the behavior of the system should not be 
biased in any way by a specific strain distribution due to 
a particular choice of the initial configuration. Further, 
the computational demands of CPMD calculations pro- 
hibit the study of loose knots because of the large system 
size required. Moreover, the initial stress must be close 
to, yet below, the breaking threshold to allow the rup- 
ture to occur within the typical time scales (< 10 ps) 
of CPMD. With these constraints in mind, we first per- 
formed classical MD simulations with the UA model. 




FIG. 5. Initial configurations of the chains containing 
N=1AA (a), 72 (b), and 50 (c) pseudoatoms, respectively. 

We chose as the initial system a chain of 144 pseu- 
doatoms, and generated the trefoil knot by adding an ap- 
propriate set of gauche defects to a linear polymer (Fig. 
panel a) . As we mentioned in section |l[ we performed 
the MD calculations by anchoring the end atoms of the 
polymer at a fixed separation distance L, which was then 
successively increased from one set of simulations to the 
next. As the knot was tightened, fewer atoms were di- 
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FIG. 6. Classical strain energy distribution along the 
chains as obtained after room temperature MD evolution and 
annealing. Shown are results for chains containing A''=50 
(solid line), 41 (dashed line), 35 (dotdashed line), and 28 (dot- 
ted line) pseudoatoms. Polymer lengths (horizontal axis) are 
scaled in order to facilitate direct comparison among different 
length molecules. 

Fig. ^ shows that decreasing the number of atoms to 
= 41 does not produce any noticeable effect; when 
A^ = 35, however, the presence of the knot begins to af- 
fect the strain distribution of the system in a significant 
way H]. As we will see later, that strain distribution 
is a general fingerprint of the presence of a trefoil; dis- 
tortions from ideal geometry tend to concentrate on the 
entrance and exit points of the knot, while its central 
portion is much less stressed. By decreasing of the num- 
ber of pseudoatoms to A^ = 28, the system becomes very 
close to the "critical" limit. The ideal chain length for 
our CPMD study thus ranges between 28 and 30 carbon 
atoms. For both of these values of A^, we performed sev- 
eral 500 ps length classical simulations at different tem- 



peratures, which were started from varying initial con- 
ditions. This was followed by 100 ps of slow simulated 
annealing to obtain a statistically meaningful set of con- 
strained equilibrium configurations for the trefoil under 
tension. We observed that, for any given length L, devi- 
ations in the final configurations obtained from different 
simulations were so small that the strain energy distribu- 
tion curves could actually be superimposed on each other 
with no noticeable difference. The initial configurations 
for the CP simulations were generated by averaging this 
set of constrained ground-state positions. 



B. CPMD calculations 




FIG. 7. Sample of initial configuration used in the ab initio 
simulations of the CisH^s molecules. The positions of the 
carbon atoms are obtained from the classical MD simulations, 
and hydrogens are added to each carbon according to the 
appropriate tetrahedral symmetry. 

The knotted minimum energy configuration obtained 
from the previous section was decorated with hydrogens 
atoms according to the ideal tetrahedral symmetry of 
carbon atoms in alkanes. Structural relaxation showed 
that this initial guess for the hydrogen positions was very 
close to those found with minimum energy geometry in 
the distorted system. In Fig. we show a sample starting 
configuration of a C2?,H^^ n-alkane. 

CP simulations were carried out with the same tech- 
nique adopted for the linear alkanes, in which external 
tensile stress was applied by constraining the end-to- 
end distance L. Results for the C28H5S and C3oHe2 



molecules were very similar. Henceforth, unless speci- 
fied, we will only refer to calculations for the N — 28 
carbon alkane. Several sets of simulations were carried 
out for L = 10.75,11.50,12.50,13.00,13.50,13.75, and 
14.00 A, respectively. For each value of L, the initial posi- 
tions were optimized (with the constraint) before heating 
the system to room temperature. We employed smaller 
timesteps and longer simulations as compared to the lin- 
ear molecules case. As explained above, this should pre- 
vent the initial conditions from biasing the behavior of 
the system. The sound velocity of polyethylene along the 
chain direction in the crystal ||l^ is about 18 km/s and 
thus an estimate of the time needed for the slowest longi- 
tudinal phonon to travel along a 30-atom chain is about 
0.25 ps; this is the lower bound for the time scale of our 
calculations. We actually performed CPMD simulations 
at least 3 or 4 times longer than this value. 
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FIG. 8. Distribution of bond lengths and angles in the 
C2SH5S alkane &t L = 13.0 A. The largest deviations from 
the equilibrium values occur at the entrance and exit points 
of the knot. 

The optimized ah initio structure is typically more fa- 
vorable by 100-120 kcal/mol as compared to the classical 
ground-state configuration. The average classical strain 
energy distributions obtained during finite temperature 
CPMD runs show a very similar shape to the classical 
ground-state distributions, but with a larger amount of 



strain located on the bonds at the entrance and exit 
points of the knot, along with a relaxation of the cen- 
tral portion of the knot. This is analogous to the results 
reported elsewhere Q and suggests that localized spikes 
in the strain distribution are necessary to allow global re- 
laxation for the remaining portions of the carbon back- 
bone. In analogy with the linear case, it is interesting to 
note that the second bond at each end is the least dis- 
torted of the whole chain. This bond is neither inside the 
knot nor the one at which the constraint force is applied. 




FIG. 9. Snapshots of CP dynamical evolution of the knot- 
ted C2SH58 alkane. For the sake of clarity we only display 
the carbon backbone. The end-to-end distance in the simula- 
tions presently shown is L = 13.5 A, and the temperature is 
T = 300 K. 

The system does not cross the dissociation barrier at 
room temperature within the time scale of the simula- 
tions for elongations smaller than L = 13.5 A. Larger 
and larger fluctuations of bond lengths, up to 30 % of 
the equilibrium values (1.54 A), are observed along the 
dynamical evolution as L increased. This is particularly 
true for the bonds located at the entrance and exit points 
of the knot; however, no bonds dissociated. In Fig^|8|we 
report an example of the C-C bond length and CCC an- 
gle distributions along the molecule for L ~ 13.0 A. The 
largest deviations from equilibrium values were about 
10 % in the lengths and 30 % in the angles; they oc- 



curred near the entrance and exit points of the knot. In 
contrast, the average deviations in the central portion of 
the knot were ca. 5 % and 18 %, respectively. 

Our data analysis indicated that monitoring the bond 
angles about the two atoms forming a bond was a 
much more effective method to determine chain rup- 
ture as compared to analysis of bond lengths. Given 
two bonded carbon atoms Ci and Cj=i+i, we refer to 
the other atoms bonded to them as Ci-i, Hi^i, Hi^2 and 
Cj+i,Hj^i,Hj^2, respectively. We define ai^n=i,?, to be 
the angles Ci-iCiHi_i, Ci-iCiHi_2 and Hi^iCiHi^2, and 
Pj,n=i,3 to be Cj+iCjHj^i, Cj+iCjHj^2 and Hj^iCjHj^2, 
accordingly. I]n'^*."=i>3 and J2n (^3,ri=i,3 equal 328.4° 
in perfectly tetrahedrical {sp^ hybridization) geometries, 
and 360° in planar (sp^) molecules. 
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FIG. 10. Upper panel: time evolution of the lengths of the 
bonds located at the entrance and exit points of a trefoil knot 
with N=28 C atoms and end-to-end distance L = 13.5 A. 
Lower panel: time evolution of the sums $ and Q of the 
bond angles about the atoms involved in the rupture (see 
text). Even when large bond-length deviations from equilib- 
rium values occurred, the chemical bond did not break unless 
a change in the hybridization status of the key atoms was 
observed. Arrows indicate that the rupture and the readjust- 
ment of the system are separated by a time interval of about 
250 fs. 

When L — 13.5 A, as shown in Fig. H, the molecule 
finally snapped during the finite temperature run at the 
exit point of the knot with the two ends then strongly 
recoiling back. In the lower panel of Fig. |l0| we report 

^22 = Z]n"22,n and 623 = Z]n/323,n aS fuUCtiouS of 

time. The fact that these quantities tended towards a 
value around 360° provides a simple, yet unequivocal, 
sign that the hybridization status of C22 and C23 atoms 



had changed from sp'^ to sp^ and, thus, that their chem- 
ical bond was broken. In the upper panel of the figure, 
we display the entrance and exit bond lengths for the 
same system. As previously discussed, very large devia- 
tions from the ground-state geometry are present for both 
bonds; however, only one of them finally ruptured. It is 
interesting to point out that, as marked by the arrows, 
the symmetric counterpart of the breaking bond needed 
about 250 fs to readjust and then oscillate around the 
equilibrium length of 1.54 A. This is in agreement with 
our estimate of the typical time scale of an alkane of this 
size. An important quantity provided by first-principles 
MD is the electron charge density of the system, which 
is particularly useful in studies focused on the evolution 
of chemical bonds. In Fig. ^ a snapshot of the charge 
density after 0.8 ps of dynamical evolution is reported. 
Between the carbon atoms C22 and C23 there is a gap of 
charge density that confirms the bond breaking. 
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break. This gave birth to three daughter radicals rather 
than two. 

The amount of strain energy that a trefoil knot stores 
in a polyethylene chain before breaking is around 300 
kcal/mol, and it has to contain at least 23 carbon atoms 
to be able distribute such load along the molecule without 
breaking. This ab initio estimate is in agreement with re- 
sults obtained via classical molecular dynamics . The 
average load per bond in the tightest possible knot is 
^ 13.3 kcal/mol, which is about 80% of the correspond- 
ing value in a linear alkane. This result, and the knot- 
induced location of the break along a chain, are in good 
agreement with known properties of macroscopic knots, 
which in turn suggests a universal behavior. 

V. CONCLUSIONS 

In conclusion, we have presented an atomic level de- 
scription of the breaking of a polyethylene-like strand, 
both with and without a knot, when subjected to ten- 
sile loading. Specifically, we find that the mechanical 
strength and the strain distribution in a polymer just be- 
fore breaking are profoundly influenced by the presence 
of a simple knot. Moreover, a knotted strand, unlike the 
unknotted one, does not break at the point where the 
tension is applied but, rather, at a point just outside the 
knot. Further, the presence of a knot weakens the rope in 
which it is tied. Our findings arc in complete accord with 
the corresponding macroscopic phenomenon exhibited by 
a rope, which suggests that some of the other properties 
of knots known to sailors and fishermen for centuries Q| 
may well have a counterpart in the nanoscopic world of 
polymers. 
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J.Hutter for fruitful discussions. Special thanks go to 
H.S. Mei for help with graphics, and to K. Bagchi and 
D. Yarne for useful suggestions. 



FIG. 11. Electronic charge density of the C2SH5S, trefoil 
at the breakpoint for L = 13.5 A and T = 300 K. The two 
atoms involved in the bond breaking are drawn with a larger 
size and a darker color. A gap in the electronic density is 
clearly observable. 

Other simulations performed using different initial con- 
ditions {i.e. heating rate, hydrogen atom optimization, 
etc.), or at larger tensile stress {L — 13.75, 14.00 A) gave 
similar results; namely, bond breaking occurs randomly 
at one of the two symmetric entrance and exit points. 
In the latter case, when L = 14.0 A, the molecule ac- 
tually broke at both points; the amount of external load 
was so high that, after the first rupture, the system was 
not able to readjust quickly enough to avoid the second 
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